Structural and dynamical heterogeneity in a glass forming liquid 
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We use the "isoconfigurational ensemble" [Phys. Rev. Lett. 93, 135701 (2004)] to analyze 
both dynamical and structural properties in simulations of a glass forming molecular liquid. We 
show that spatially correlated clusters of low potential energy molecules are observable on the time 
scale of structural relaxation, despite the absence of spatial correlations of potential energy in the 
instantaneous structure of the system. We find that these structural heterogeneities correlate with 
dynamical heterogeneities in the form of clusters of low molecular mobility. 
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Over the last decade, the identification and study of 
dynamic heterogeneity (DH), especially in computer sim- 
ulations, has added an important new dimension to our 
understanding of complex relaxation in glass forming liq- 
uids P, 0- DH refers to the emergence and growth 
of spatially correlated domains of mobile and immobile 
molecules as temperature T approaches the glass transi- 
tion temperature T g . A question that dominates research 
on DH concerns its connection to the structure of the 
liquid: What local configurational properties influence 
whether a given molecule is mobile or immobile? 

Recent work by Widmer-Cooper, Harrowell and 
Fynewever [j| has shown conclusively that a structure- 
dynamics connection must exist at the molecular level. 
To do so, they use an "isoconfigurational (IC) ensem- 
ble" 0, \3 13 , a set of microcanonical molecular dynam- 
ics (MD) trajectories in which each run starts from the 
same initial equilibrium configuration, but with molecu- 
lar momenta chosen randomly accordingly to a Maxwell- 
Boltzmann distribution. The result is a set of trajectories 
lying on the same energy surface, and evolving away from 
their common initial point in configuration space. They 
then define and evaluate the "dynamic propensity" : the 
average, in the IC ensemble, of the squared displacement 
of a molecule at a time comparable to the structural re- 
laxation time. They show that DH is observed in this 
approach, in the form of increasing spatial correlations 
of the dynamic propensities in a glass forming liquid as 
T — > T g . Since the influence of the initial momenta is 
averaged over, the observed spatial correlations must be 
configurational in origin. 

The strength of Ref. Q is that it exposed the fea- 
tures of DH that are structural in origin, without needing 
to determine what structural properties are responsible. 
Other studies have worked towards explicitly identifying 
structural correlators to dynamics. A number of works 
over the past decades have identified relationships be- 
tween average structural properties (especially free vol- 
ume and measures of symmetry in the local molecular en- 
vironment) and bulk dynamics; see e.g. Refs. 0,011)0- 
More recently, several studies have sought a correlation at 
the microscopic level, e.g. between local free volume and 
local mobility, with more success in some systems 0, |nj 
than in others A notable absence of correlation be- 



tween the local volume and the local Debye- Waller factor 
has been reported recently Insights have also been 
realized using local measures of symmetry to elucidate lo- 
cal mobility, e.g. icosohedral order [itLI14| : or the recent 
work of Shintani and Tanaka in which less mobile regions 
were shown to correlate well to structured, crystal-like 
domains 0] . Recently, the local Debye- Waller factor has 
been shown to correlate to the dynamic propensity 
firmly establishing the connection between local dynam- 
ics at short and long times. 

Notably absent from the list of local structural quanti- 
ties that correlate well to local dynamics is the poten- 
tial energy. It has been shown that a molecule with 
a low potential energy will be less mobile, on average, 
than one with a high potential energy 0]. However, the 
variance around this average trend is comparable to or 
larger than the trend itself, making it impossible to pre- 
dict what a given molecule will do based on its potential 
energy. Careful time averaging and the use of in- 
herent structures 0] has been shown to yield little gain 
in correlation. This is both perplexing and disappoint- 
ing. Perplexing since we know that the average potential 
energy strongly influences the average dynamics of the 
system . And disappointing because the potential en- 
ergy is a natural observable to correlate to dynamics, es- 
pecially given the interest in the analysis of glass-forming 
systems using the potential energy landscape |l9j . 

In this Letter, by expanding the application of the IC 
ensemble to structural quantities, we show that it is pos- 
sible to identify heterogeneities of the potential energy 
that correlate well to dynamic heterogeneities, in a liq- 
uid where no useful correlation is discernible from the 
instantaneous properties of the system. Our work shows 
that emergent "structural heterogeneities" , that develop 
in tandem with dynamic heterogeneities, exist and can 
be observed in a glass forming liquid. 

Our results are based on molecular dynamics simula- 
tions of N = 1728 water molecules interacting via the 
ST2 pair potential 20] . Much is known about this simu- 
lation model, both in terms of thermodynamic and trans- 
port behavior, making it a good candidate for a detailed 
analysis of structure-dynamics relationships in a model 
3D molecular liquid |1 IH |H . We study three T 
(350, 290 and 270 K) all at density p = 0.83 g/cm 3 . At 
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FIG. 1: (a) Kt and (b) Cv as a function of T along the p = 
0.83 g/cm 3 isochore. Data are derived from the simulations 
described in Ref. 12311. 



this p, the hydrogen bond network in this water model 
is more prominent than at other p, suggesting that lo- 
cal energetics may have a particularly strong influence 
on dynamics. This p is also a convenient choice because 
the isothermal compressibility Kt (which is related to 
the strength of static density fluctuations) is decreasing 
with T (Fig.QJ. This ensures that any DH that emerges 
as T decreases will not be due to the growth of conven- 
tional density fluctuations. At the same time, we note 
that the isochoric specific heat Cy (and hence the mag- 
nitude of fluctuations of the system energy) is increasing 
to a maximum, so spatial variations of potential energy 
may be occurring. These conditions therefore provide a 
promising regime in which to seek connections between 
energy and dynamics at the molecular level. 

In all our simulations, molecular interactions are cut 
off at a distance of 0.78 nm, and effects of those at longer 
range are approximated via the reaction field method. 
All simulations have constant N and volume V, and use 
a 1 fs time step. At each T, we first conduct a stan- 
dard MD run where T is maintained near the desired 
value using the method of Berendsen, et al j2f|. To en- 
sure equilibration, these runs are carried out for twice the 
time required for the mean squared displacement to reach 
1 nm 2 . We then use the configuration at the end of the 
equilibration run as the initial configuration for generat- 
ing M = 1000 runs of an IC ensemble. Note that both 
the linear and angular momentum of each molecule is 
randomized, using the Maxwell-Boltzmann distribution 
appropriate to the T of the initial configuration. Each 
run is carried out in the microcanonical ensemble for a 
time t = 308 ps (350 and 290 K) or 937 ps (270 K). 

Let r 2 (i,k,t) be the squared displacement of the 
O atom of molecule i at time t in run k of 
an IC ensemble. The system-averaged and IC- 
ensemble- averaged mean squared displacement, (r 2 ) — 
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rameter, a = [(3(r 4 ))/(5(r 2 ) 2 )] 



and non-Gaussian pa- 
— 1 are shown in Fig. [21 
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FIG. 2: (a) Mean squared displacement (r 2 ), (b) non- 
Gaussian parameter a, (c) mean cluster size for least mobile 
molecules, SV, and (d) mean cluster size for tightly bound 
molecules, S u , all as a function of t. From left to right in 
(a), the curves are for T = 350, 290 and 270 K. In (b-d), the 
curve with the left-most maximum is T — 350 K, the middle 
maximum is 290 K, and the right-most is 270 K. 



Both show the characteristic pattern of a glass forming 
liquid in which DH occurs, (r 2 ) develops a plateau at low 
T indicating the onset of molecular caging, and a displays 
an increasingly prominent maximum as T decreases. We 
have checked that the curves in Figs.Efa-b) are, at large 
t, within error of those found in a conventional approach, 
where an average is taken over multiple time origins of 
a single, long MD run. This suggests that the particular 
initial configurations chosen are at least approximately 
representative of the equilibrium behavior. 

The dynamic propensity of each molecule is the value 

of (r 2 )ic = A/ -1 Y^k=i r2 (*> t) when t is on the or- 
der of the structural relaxation time It measures 
the propensity for molecule i to undergo a given dis- 
placement, given its starting point in the initial config- 
uration, rather than indicating what the molecule will 
do in any particular run. Here, we extend the use of 
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FIG. 3: Dynamical heterogeneity (left panels) and structural 
heterogeneity (right panels) in the initial configuration at T = 
350 K (top panels), 290 K (middle panels) and 270 K (bottom 
panels). To make each panel, the values of (rf)i C (or («i)i c ), 
evaluated at the time of the maximum of S u , are assigned 
to each molecule in the initial configuration. These values 
are sorted and assigned an integer rank Ri from 1 to N, from 
smallest to largest. Each O atom is then plotted as a sphere of 
radius i? min exp{[(Ri - N)/(l - N)] log(_R max /7? min )}, where 
Rmmx = 0.14 nm and i? m i n = 0.004 nm. The result repre- 
sents the rank of {r?)i c or (ui)i c on an exponential scale, such 
that the largest spheres on the left represent the least mo- 
bile molecules, and the largest on the right the most tightly 
bound. Note that hydrogen atoms are not shown. 



the IC ensemble to study structural properties as well 
as dynamics. We focus on u(i, k, t), the contribution of 
molecule i to the instantaneous potential energy of the 
system at time t in run k of an IC ensemble. Specifi- 
cally, Ui = Ylj=i wn ere <pjj is the pair potential en- 
ergy of molecules i and j. Analogous to (rf)i c we define 
(wj)ic = M^ 1 '^2 k _ 1 u(i,k,t). At a fixed t, (u,)i c mea- 
sures the propensity of molecule i to have a given value 
of the potential energy, given its starting point in the 
initial configuration. 

First we test for the occurrence of DH, by evaluating 
(rf)i c for each molecule as a function of t, and examin- 
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FIG. 4: (a) u° and (b) (wi)ic, both versus (rf )i c at each T. 
The times indicated, at which (ui)i c and (rj)i c are evaluated, 
correspond to the maxima of S u - To ease comparison, in (a) 
the cloud for 290 K has been shifted upward by 50 kj/mol, 
and that for 350 K by 130 kj/mol; in (b) the cloud for 350 K 
has been shifted down by 4 kj/mol. 



ing the s pat ial arrangement of this quantity via a cluster 
analysis |26l |2?J • For reasons that will become clear be- 
low, we focus on the least mobile molecules, specifically 
the subset having the lowest 10% of (rf )j c values. Clus- 
ters are defined by the rule that two molecules of this 
subset that are also within 0.35 nm of one another (the 
position of the first minimum of the 0-0 radial distri- 
bution function) in the initial configuration are assigned 
to the same cluster. The number-averaged mean cluster 
size S r is then found from (1/N C ) ^ n nC(n), where C(n) 
is the number of clusters of size n, and N c is the total 
number of clusters. Fig. |2Jc) shows the t dependence of 
S r at each T. At small t, S r has a value consistent with 
a random choice of 10% of the molecules, approximately 
S r = 1.27, indicating no spatial correlations in (rf)i c . 
However, at intermediate times corresponding to the on- 
set of structural relaxation, a maximum occurs, indicat- 
ing significant clustering of the least mobile molecules. 
At large t, S r begins its return to the uncorrelated value, 
and the DH dissolves. The morphology of the DH ob- 
served near the maxima in Fig. Hfc) is illustrated in the 
left panels of Fig. [21 The correlated domains of larger 
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spheres in Fig. [3| indicate the locations of the large clus- 
ters that generate the maxima in Fig. |5fc) . 

We then carry out exactly the same analysis, but us- 
ing the lowest 10% of (itj)i C values; this selects the subset 
of molecules with the greatest propensity to be "tightly 
bound." The mean cluster size for this subset, S u , shows 
a similar behavior to S r [Fig. |2td)]. Note that in the 
limit t — > 0, we have (iti); c u°, where u° is the instan- 
taneous potential energy of each molecule in the initial 
configuration. The behavior of S u at small t confirms 
that u° shows essentially no spatial correlation at any T. 
And yet, on the same time scale as the signature of DH 
is observed in S r , a maximum also occurs in S u . The 
spatial correlations of molecular potential energies (i.e. 
structural heterogeneities) that generate the maxima in 
S u are illustrated in Fig. [31 The morphology of the two 
types of heterogeneity shown in Fig. [3] is strikingly sim- 
ilar: there is a clear correlation between regions with a 
propensity to be the least mobile, and a propensity to be 
tightly bound. It is worth bearing in mind that these two 
types of heterogeneity are defined independently. Neither 
the time scale on which the structural heterogeneity oc- 
curs, nor the definition of the structural clusters depends 
on dynamical information. 

The absence of a useful correlation between u° and the 
dynamic propensity is demonstrated by a scatter plot of 
u° against (ff)i c , evaluated at the time of the maximum 
of S u at each T [Fig.^Ja)]. The best correlation is found 
at high T, but even here some molecules in the lowest 
10% of u° have values of (rf )i c comparable or even larger 
than the mean of (r?)j c - At lower T the correlation only 
gets worse: molecules with the lowest values of u° are 
found across the entire spectrum of (rf ) ic values. We 
have also tested if the correlation improves if we replace 
u° in Fig.^Ja) with its value in the inherent structure of 
the initial configuration, but it does not. 

Fig- Elk) shows a scatter plot of (iti)j c versus (r|)i c at 
the time of the maximum of S u at each T. Here we see 
the reason why we have focussed on the least mobile and 
most tightly bound molecules: it is at the lower end of 
these scatter plots that the points are most easily dis- 
tinguished from the overall population. In contrast, the 
correlation between the most mobile and least tightly 
bound molecules is little improved over that in Fig.|3Ja). 
We have also examined the clusters formed by the 10% 
most mobile and least tightly bound subsets of molecules. 
The most mobile molecules also form clusters at interme- 
diate times, although the strength of the effect is weaker 
than for the least mobile molecules. Interestingly, the 
least tightly bound subset shows a decreasing tendency 
to cluster as T decreases, with the maximum in S u be- 



coming difficult to discern at the lowest T. Whether 
this is a generic behavior, or a feature of this particular 
network-forming liquid requires further research, which 
is currently underway. One possibility is that molecular 
hopping is emerging as an important mode of transport 
at low T, which might weaken the energy-dynamics cor- 
relation for the most mobile molecules [23 |. 

We can therefore make the following statements about 
the relationship of local potential energy and local dy- 
namics in this system: Despite the absence of a corre- 
lation between the instantaneous potential energy of a 
molecule and its subsequent displacement, a molecule 
that has a propensity to be tightly bound (as evalu- 
ated within the IC ensemble) also has a propensity to be 
among the least mobile. On the other hand, a propen- 
sity to be mobile and a propensity to be loosely bound 
do not correlate at low T, indicating that other factors 
are controlling dynamics at this end of the spectrum. 
The IC ensemble thus makes possible a dissection of the 
energy-dynamics relationship, under conditions where no 
conclusions can be drawn from instantaneous structural 
information. The approach also exposes the morphol- 
ogy of the structural heterogeneity present in a single 
instantaneous configuration. As Ref. Q first showed, av- 
eraging over the IC ensemble therefore acts as a powerful 
amplifier of the subtle configurational signal that must 
exist instantaneously, but which is overwhelmed by the 
"noise" imposed by molecular momenta. While we have 
only addressed the behavior of the potential energy in 
this work, we expect that valuable information concern- 
ing many supercooled liquids can be extracted from a 
wide range of configurational quantities (e.g local vol- 
ume, local symmetry) when evaluated at the local level 
in the IC ensemble. 

Our work also demonstrates that the emergence and 
growth of heterogeneity in a glass forming liquid is not 
solely a dynamical phenomena, but also structural, and 
that the two types of heterogeneity can be analyzed 
within the same framework by using the IC ensemble. 
In particular, we note that the size of the structural and 
dynamical heterogeneities observed here are approaching 
the system size at our lowest T. It may be that the ability 
of the IC ensemble to discern subtle spatial correlations 
will help clarify the nature of the growing length scales 
of dynamic and structural heterogeneity as T — > T g . 
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